function y = eval_Kij(xi, xj)
    y = 1/(xj - xi).*quad(@obj, 0, 1);
    function y = obj(x)
        y = phi0(x).*phi1(x);
    end
end